Kelvin-Helmholtz instability in high-energy heavy-ion collisions 
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The dynamical development of collective flow is studied in a (3+l)D fluid dynamical model, with 
globally symmetric, peripheral initial conditions, which take into account the shear flow caused by 
the forward motion on the projectile side and the backward motion on the target side. While at 
y^ijvjv = 2.76^4 TeV semi-peripheral Pb+Pb collisions the earlier predicted rotation effect is visible, 
at more peripheral collisions, with high resolution and low numerical viscosity the initial development 
of a Kelvin-Helmholtz instability is observed, which alters the flow pattern considerably. This effect 
provides a precision tool for studying the low viscosity of Quark-gluon Plasma. 

PACS numbers: 24.85.+p, 24.60.Ky, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 

Global collective observables are becoming the most es- 
sential in ultra-relativistic heavy ion reactions When 
we want to extract information from experiments, both 
on the equation of state (EoS) and the transport prop- 
erties of matter 0, @1 , we have to invoke a realistic de- 
scription with a fully (3+l)D dynamical evolution at all 
stages of the reaction, including the initial state. 

It is important to note that the phase transition to 
quark-gluon plasma (QGP) and consequent fluctuations 
may enhance the collective behavior of the system 
For the fluid dynamical (FD) initial state we must have 
a system that is close to local equilibrium; thus, at high 
energies the transition to QGP has to happen earlier than 
the formation of the locally equilibrated initial state. 

The (3+l)D, relativistic FD model we use to describe 
energetic heavy ion reactions is well established and de- 
scribes the measured collective flow reliably p| @ . We 
use the Particle in Cell (PIC) method in which an Eu- 
lerian grid contains a very large number of Lagrangian 
marker particles that move with the matter. This method 
enables us to follow the motion of the fluid with good 
precision. At Large Hadron Collider (LHC) energies in 
these calculations we observed a significant rotation of 
the QGP fluid in peripheral collisions, which leads to ob- 
servable consequences @. 

Our detailed studies indicate the development of an in- 
teresting phenomenon, namely the beginning of a physi- 
cal instability. In peripheral collisions, in the transverse, 
[y, z] plane a non-sinusoidal instability starts to develop. 
We can visualize this by coloring the markers in the pro- 
jectile (upper) side blue and the target (lower) side red. 
Initially the dividing surface is a plane. As time proceeds, 
the markers (which indicate the location of conserved 
baryon charge) move, and the dividing surface becomes 
a wave, which resembles the start of a Kelvin-Helmholtz 
(KH) instability as shown in Fig. [T]in the [x, z], reaction 
plane, i.e., \y\ < 1 cells (\y\ < 0.7 fm). 



Initially, at 1.5 and 3.0 fm/c, we can see two shorter 
wavelength KH instabilities, which then dissolve and are 
fed into a longer wavelength instability. 

The non-sinusoidal behaviour of the instability at later 
times is not obvious as the development of instability and 
the spherical expansion compete with each other. 

The density of the central zone decreases rapidly so 
that the matter "freezes out" and the fluid dynamical 
description breaks down. As a consequence we do not 
see the vorticity sheet rolling up as in a fully fledged KH 
instability. In more central collisions the dividing surface 
nearly remains a plane. The usual reason for the KH 
instability is a "shear flow" , where in a fluid layer there is 
a large velocity gradient. Thus, the origin and the energy 
source of this phenomenon is in the initial configuration 
and the initial velocity distribution, which are correctly 
represented in our initial state model Q- 

In our computational fluid dynamics (CFD) calcula- 
tion the initial state model -based on longitudinally ex- 
panding flux tubes or streaks is used @. In non- 
central collisions only a part of the original nuclei inter- 
act. These are in the participant zone where the streaks 
develop. Spectator nucleons on the two sides are not 
participating in the reaction. 

The participant streaks are formed by the color charges 
arising from the projectile and target nuclei after these 
have penetrated through each other. The chromo-electric 
field, characterized by the string tension, slows down the 
expansion of the ends of the streaks. Our FD initial state 
is a configuration where the matter is stopped within 
each streak, while streaks expand independently of each 
other. Thus, this model is applicable streak by streak 
and the momentum of the streaks varies, especially for 
the peripheral streaks where the asymmetry between the 
projectile and target contribution to the participant mat- 
ter is the biggest. So, the streaks at the projectile and 
target sides move in the beam (z) direction with sub- 
stantial velocity difference. This generates the shear flow 
configuration. 
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FIG. 1: (color online) Growth of the initial stage of Kelvin- 
Hehnholtz instability in a 1.38A + 1.38A TeV peripheral, 
b — 0.76 max , Pb+Pb collision in a relativistic CFD simu- 
lation using the PIC-method. We see the positions of the 
marker particles (Lagrangian markers with fixed baryon num- 
ber content) in the reaction plane. The calculation cells are 
dx — dy — dz = 0.4375fm and the time-step is 0.04233 fm/c 
The number of randomly placed marker particles in each fluid 
cell is 8 3 . The axis-labels indicate the cell numbers in the x 
and z (beam) direction. The initial development of a KH 
type instability is visible from t — 1.5 up to t = 7.41 fm/c 
corresponding from 35 to 175 calculation time steps). 

The aspect ratio of matter in the reaction plane ([x, z\- 
plane) becomes more elongated with increasing impact 
parameter, 6, as the height of the participant profile, L, 
becomes smaller, L — (2R — 6), for nuclei of radius R, 
and the streaks are becoming longer due to the smaller 
effective string tension 7]. Thus, the aspect ratio for 
b = 0.5& max (where 6 max = 2R) is [1 : 1.5], while for 
b = 0.7& max it is [1 : 3]. Of course this aspect ratio de- 
pends on the initial state model, and some of these do 
not take into account the longitudinal expansion before 
thermalization, and even less the dependence of the ex- 
pansion on the effective string tension. 

In a heavy ion reaction the projectile edge of the par- 
ticipant domain moves almost with the velocity of the 
projectile, u, while the target side moves with the target 
velocity, —u. At high energy this difference provides con- 
siderable shear in the velocity fields. At the same time in 
the initial state model [7| the initial transverse velocity 



is zero for all fluid elements. A low [x : z]-profile makes 
it possible to develop a typical shear-flow configuration 
and thus, there may be a possibility to form the initial 
stages of a KH instability. 



II. PHYSICAL CONSIDERATIONS 

A. Growth of the KH instability 

The growth of a small initial KH instability in an ideal- 
ized shear-now configuration can be described in a rather 
simple way. From ref. [|| sect. 3 it follows that the shear 
flow starting from a small sinusoidal perturbation in in- 
compressible and inviscid now, the perturbation will grow 
exponentially, oc exp(si), (Q 3.15), where t is the time 
and s is proportional to the wave number k, 



kV 



(Q 3.28) where ±V is the characteristic velocity of the 
upper/lower sheets which is somewhat less than the pro- 
jectile or target velocity. 

Thus, the largest k or shortest wavelength will grow 
fastest. Also, increasing the beam energy (i.e. increasing 
V) will also lead to increased development of turbulence! 

For a Pb+Pb reaction, R = 7 fm, b max = 14 fm, we 
study the impact parameters b = 0.5 or 0.76 max . For 
these collisions the typical transverse size of the initial 
shear flow is 

L = (2R — b) = 7.0 or 4.2 fm . 

The typical calculation cell size is dx = dy = dz — 0.35 
fm. The beam directed, longitudinal length of the initial 

state is 

l z = 10.5 or 13.1 fm , (1) 
and the minimal wave number is 

k = 2n/l z = 0.598 or 0.479 fm -1 

respectively. 

For the scaling analysis of instabilities we need the di- 
mensionless numbers constructed from the typical length, 
L, and speed, V. The Reynolds number is Re — 
VL/u, where v is the kinematic viscosity. So, for a 
peripheral heavy ion collision with impact parameter 
b, the characteristic length is L and V is the veloc- 
ity of the top/bottom layer, V = \u\. In exactly cen- 
tral collisions u = 0, while with increasing impact pa- 
rameter, u = ±0. 26, 0.34, 0.36, 0.39, 0.43, 0.42, 0.39 c for 
b = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 & max , respectively. No- 
tice that due to the geometry and the minimal string 
diameter of 1 fm, the increase of this velocity does not 
reach the beam velocity, so for typical peripheral colli- 
sions V = \u\ ~ 0.4c. 

In the simplest incompressible and inviscid flow ap- 
proximation the amplitude of of the starting turbulence 
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would double in 2.90 or 3.62 fm/c for b = 0.5 or 0.76 max . 
The growth of instability is very fast in this approxima- 
tion and it increases with the beam energy (beam ve- 
locity) and with the wave number. The typical reaction 
time in a heavy ion collision exceeds the time needed to 
double the amplitude of an initial instability. 

At the same time we also observe that at smaller im- 
pact parameters the development of an instability is not 
seen in our calculations, see Fig. [2j Thus, we have to 
conclude that the role of viscosity is decisive as a large 
viscosity will decrease this growth rate and may eliminate 
the possibility of the KH instability. 

B. Formation of critical length KH instability 

While perturbations with larger wave number (shorter 
wavelengths) may grow faster, there is a critical minimal 
wavelength beyond which the perturbation is stable and 
able to grow. Smaller wavelength perturbations tend to 
decay into random thermal fluctuations. This situation 
is analogous to the phase transition dynamics via homo- 
geneous nucleation where the formation of critical size 
bubbles or droplets is required to start the phase transi- 
tion 0. 

This aspect of turbulence formation was first discussed 
by Kolmogorov [9] for flow in the "inertial range" where 
the effects of viscosity are still negligible. These consid- 
erations are applicable until the viscosity does not have a 
significant effect on the formation of vortices. The min- 
imal stable wavelength for a starting instability is given 
by the Kolmogorov length scale, which is the smallest 
scale of turbulence. Dominant and increasing viscosity 
results in increasing critical vortex size, AkoI- Smaller 
perturbations are unstable. 

The average rate of energy dissipation, e, per unit mass 
and unit time, is associated with the decay of an eddy of 
size I and characteristic speed, vi into two smaller ones, 
in time t; = l/vi. It follows then that e ~ vf/ti — vf/l. 
The decay (or formation) time of the small size eddy, 
ti, can be compared to the viscous diffusion time of a 
perturbation of size I, which is tf ts - = I 2 /v. Equating the 
two estimated characteristic times provides the minimal, 
Kolmogorov length, \ko\- 

A Ko1 = [v 3 /e]V 4 . (2) 

Here v is the kinematic viscosity of the fluid, v = rj/ p — 
r]/(nmB), where r\ is the shear viscosity, p is the mass 
density, n is the baryon charge density and ms is the 
characteristic mass falling on unit net baryon charge in 
QGP. 

Even if the average rate of energy dissipation, e, is pro- 
portional to the viscosity, this dependence is linear, so 
the critical vortex size is still increasing with increasing 
viscosity. The key question is: can a critical size vor- 
tex be formed in a heavy ion collision? Lower viscosity 
and higher energy or energy dissipation may enable the 
formation of a critical size or larger vortex. 



Kolmogorov's theory also provides an energy distribu- 
tion spectrum for small vortices or whirls in nearly per- 
fect fluids. The energy density spectrum in terms of the 
scale of the vortices, A, is proportional to A -1 / 3 , i.e., it 
is lower for larger vortices. Large vortices may generate 
smaller ones, until we reach the viscous limit, AkoI, where 
the vortices are becoming just thermal fluctuations. 

In non-relativistic flow the mass flow is identical with 
the flow of massive particles while the flow energy and 
the thermal energy are negligible compared to the rest 
mass of the fluid. At ultra-relativistic energies this is not 
the case and the separation of flow (inertial) energy and 
the random thermal energy is not a trivial question. In 
order to follow the classical concepts turbulence and its 
development we will follow Eckart's definition of flow ve- 
locity, where flow is bound to the conserved net baryon 
charge and we will assume that the flow of the average of 
all quarks can be characterized by this velocity. The ob- 
served constituent quark number scaling of collective flow 
observables for different hadrons supports this approach. 

We focus on the description of the initial stages of the 
development of turbulent instability in the central zones 
of the collision at a period just after the formation of the 
locally equilibrated FD initial state. We will estimate 
the corresponding collective mass of the partonic mat- 
ter in QGP per unit baryon charge, ms, at this stage of 
the reaction. Before the collision at the LHC each nu- 
cleon has 1.38 TeV energy. Local equilibration is reached 
when, for b = 0.76 ma x impact parameter, in the first 3 
fm/c time after initial equilibration the average temper- 
ature is T w 400 -7- 600 MeV, the average entropy den- 
sity is s ~ 150 -7- 440 fm -3 , the average baryon charge 
density isri = 0.1-^0.16 fm~ 3 and the average internal 
energy is T(s/n) = Ta = e/n w 1.1 ~ 1.3 TeV/nucl., 
where a is the specific entropy per unit baryon charge. 
At the initial stages of QGP flow the remaining energy is 
shared between collective flow energy and particle mass 
of all constituents, plus a smaller amount may go to pre- 
equilibrium emission of high energy particles. Thus, the 
effective mass can be estimated as mg « 100 GeV per 
unit baryon charge. 

By late stages of the reaction the kinetic energy of 
the flow increases substantially, while the dissi pat ion in- 
creases the thermal energy by about 4 6% [10(. The 
total hadron multiplicity increases by about an order of 
magnitude, so the effective mass per net baryon charge 
is of the order of tub ~ 10 GeV. We are interested in 
the case of initial QGP in local equilibrium when the KH 
instability could start. 

The shear viscosity is temperature dependent with a 
sharp minimum at the critical point of the phase transi- 
tion between hadronic matter and QGP [3], and reaching 
unity at the initial hot, compressed QGP at about 4T C : 

s 

which for the initial QGP gives 

7] = sh « 30 ^ 158 GeV/(fm 2 c), 
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and about 10 times less for the minimal viscosity. Con- 
sidering that the initial effective mass density may vary 
between, p = nms ~ 10 -f- 16 GeV/(fm 3 c 2 ), the corre- 
sponding kinematic viscosity is 



flow energy to heat is 



= 2.5 ~r~ 16 fm< 



and about 10 times less for the minimal viscosity. 

The corresponding Reynolds number is Re = 0.3 — 1 
(for "rj/s = 1"), and if we choose the minimal viscosity 
(V s = 0- 1 "), then Re = 3-10. These are small Re val- 
ues for turbulent flow in general, but the KH instability 
can also appear for small Re 8] . 
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FIG. 2: (color online) Time evolution of the flow in 1.38T + 
1.38A TeV peripheral, 6 = 0.56 m ax Pb+Pb collision. The 
calculation cells are dx = dy = dz = 0.585 fm and the time- 
step is 0.08466 fm/c. The number of randomly placed marker 
particles in each fluid cell is 8 3 . In contrast to Fig. [1] the 
KH instability does not develop during the initial ~ 8 fm/c 
time, due to the increased numerical viscosity and the similar 
length and height of the initial state. Thus the sides of the 
fluid provide a stiffer formation and the system rather rotates 
as a solid body, instead of forming an expanding turbulent, 
rotating shell. Due to angular momentum conservation the 
rotation slows down as the system expands. 

If we have a perfect fluid the flow is adiabatic, there is 
no dissipation, so AkoI ~ (0/0) ! The specific dissipated 



e = e/p oc To I p oc v, 



(3) 



where the dot indicates the proper-time derivative, e = 
d t e is the change of energy density with time, T is the 
temperature and a is the proper time derivative of the 
specific entropy density, a. Thus if v — > then Ak i — > 
0. So, the minimal size will grow from zero if the viscosity 
grows. With finite viscosity one can have a large minimal 
size, so that the turbulence can not develop within the 
given length of the system. As in the final expanding 
stages of the QGP fluid the viscosity increases [3(, the 
minimal eddy size, AkoI, will be larger, so initial smaller 
length instabilities will disappear. 

A good example for the formation of a minimal size 
eddy can be observed in a two-component Fermi-gas (e.g. 
Li-6), which forms a super-fluid at low temperatures in a 
rotating magnetic trap [111 ]. If we reach a limiting rota- 
tion frequency, a small eddy may develop in the central 
region of the cylindrical trap when the energy of a small 
critical size eddy will be sufficient to balance the vis- 
cous dissipation, and this eddy may become stable. One 
also needs a given viscosity or scattering length to form 
this eddy in the middle, thus the eddy first appears at a 
given finite scattering length! See Fig. 14 of ref. [ll|a. 
The minimal KH instability has similar minimal size be- 
haviour, although the shear-flow geometry is different; 
one needs a minimal torque from the boundary condi- 
tion and a minimal viscosity to form a critical eddy. 

Let us make a very simplified estimate for the size of 
the smallest possible eddy in a heavy ion collision. For 
finite shear viscosity the energy dissipation per unit mass 
and unit time, e = e/p, depends on the viscosity as well as 
the flow pattern. The characteristic shear V/L depends 
on the impact parameter, so that (V/L) 2 = 0.0038 
0.0086 (c/fm) 2 for b = 0.5^0.76 max , respectively. For the 
ideal shear flow the dissipated energy with the minimal 
viscosity ("77/s = 0.1") is (see [l2j sect. 16 or [l3[ sect. 
1.6): 



V 

T 

11 26 MeVc/fm 4 for b = 0.5b lr 

61 138 MeVc/fm 4 for b = 0.76 n 



(4) 



Then the energy dissipation, with an estimated average 
unit mass density of p = 13 GcV/ (fm 3 c 2 ), gives the rate 
of energy dissipation e = e/p and thus the Kolmogorov 
length is 



AkoI — 



2.1 
1.4 



5.4 fm for b - 
3.6 fm for b ■■ 



0.56 
0.76 max 



(5) 



In peripheral heavy ion collisions the KH instability can 
develop only if 



L > A 



Kol 



(6) 
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Thus, comparing the above values of the Kolmogorov 
length scale to the length of the initial state in the beam 
direction, eq. (JT]), we can see that at b = 0.76 max we may 
have a possibility to initiate a KH type instability in a 
heavy ion collision, while at smaller impact parameters 
this possibility is marginal. This is enhanced by the fact 
that viscosity and the Kolmogorov length increases with 
expansion, so the time-slot where the KH instability may 
develop is reduced for more central collisions. 

Based on eqs. (|2|4|) we see that in our situation the 
Kolmogorov length is proportional with the square of 
viscosity, AkoI oc rj 2 , so if the viscosity doubles, the Kol- 
mogorov length, eq. ([5]), increases by a factor of four, so 
it can reach 20 fm, which exceeds the longitudinal system 
size, and hinders the development of the KH instability 
up to about 8 fm/c. This is illustrated in Fig. [2] with 
the increased numerical viscosity. The change of marker 
particle resolution does not influence the disappearance 
of the KH instability. 

The characteristic geometry of the KH instability (ap- 
proximately two planes close to each other) require that 
L < R or b > _R, and so b may vary in the interval be- 
tween 0.5& max and 0.86 max where 6 max = 2R. At larger 
impact parameters L becomes too small and the general 
applicability of a continuum approach becomes question- 
able. So, this parameter cannot vary too much, and con- 
dition (j6]) requires to have a viscosity smaller than a lim- 
iting value, 

v c ~ 5 fmc, 

which is satisfied by the low viscosity QGP, even at higher 
than critical temperatures. 



III. FLUID DYNAMICAL MODEL 
PREDICTIONS 

We have performed CFD simulations with the PIC so- 
lution method [f| H| where the equations of relativistic 
FD were solved for a perfect quark-gluon fluid. At the 
same time the numerical method, due to the finite grid 
resolution, led to dissipation and thus to entropy pro- 
duction, which has been analysed 10] . From the entropy 
production we could determine the corresponding "nu- 
merical viscosity" , and this was approximately the same 
as the estimated, low viscosity of the quark gluon plasma. 
To avoid double counting, i.e., over counting of viscous 
dissipation, we did not add additional viscous terms to 
our CFD model simulations 1 . 



Our numerical model predictions confirm the previ- 
ously presented physical conclusions, and these show a 
developing instability which is visible at b = 0.7& max 
(Fig. [T]), but at the smaller impact parameter the KH 
instability is weak and does not change with increasing 
resolution, see Fig. [3] Increasing the resolution by 67% 
at b = 0.76 max increases the amplitude of the KH insta- 
bility wave by 57%. The final KH instability amplitude 
(1.1 fm) is 6% of the final profile height (and 16% of the 
initial one). 

In the case of heavy ion collisions the special geometry, 
i.e., the shape of the participant zone, is also hindering 
the development of the instability because the more ex- 
tended side-walls at smaller impact parameters have a 
stabilising effect against the KH instability. 

We can observe in our calculations that the initial si- 
nusoidal wave shape will become asymmetric in standard 
plane shear flow (see 3 Fig. 3.3 or [3] Fig. 2), especially 
at points of accumulating vorticity. 

The PIC method has a particular advantage in study- 
ing the KH instability. The numerical viscosity is set by 
choosing a calculation grid size, which provides the es- 
timated small dissipation of the viscous QGP fluid. At 
the same time the PIC method has a large number of 
marker particles in each Eulerian fluid cell. Their num- 
ber can vary, and can be orders of magnitude higher then 
the number of Eulerian fluid cells forming the fixed cal- 
culation grid. Thus the motion of the marker-particles 
provides a fine resolution and can follow the dynamics 
of the flow more accurately. In this method the marker 
particles provide an accurate tracing of the initial devel- 
opment of the KH instability. 

In connection with our CFD solution one has to men- 
tion that the numerical viscosity of our model calculation 
is small, r//s = 0.1, based on the small Eulerian cells. In 
addition the number of initial marker particles per nor- 
mal density cell was changed from 3 3 = 27 to 9 3 = 729 so 
that the higher Lagrangian resolution allows for a more 
accurate description of the instability. 

In case of heavy ion reactions the flow is not stationary 
and the shear flow geometry is only present in the initial 
state. Later, due to the large pressure of QGP the plasma 
explodes and expands radially in a way that the final flow 
pattern at freeze out is close to spherical albeit somewhat 
elongated longitudinally and in the reaction plane (in- 
direction) due to the dominant elliptic flow. 



In numerical, CFD simulations of instability, it is important 
to study the dissipativc effects of both the numerical viscosity 
and the physical viscosity. A finite grid resolution in the CFD 
solutions leads to the absorption of the shorter wavelength and 
high frequency fluctuations, and the energy of these fluctuations 
is converted into heat. There exist solution methods, with finite 



computational grid resolution, which enforce entropy conserva- 
tion for perfect fluids. These solutions then appear to be perfect 
adiabatic fluid flow solutions. However, this is misleading, such a 
numerical solution still absorbs high frequency small wavelength 
perturbations, while the energy of these is converted into large 
wavelength fluctuations. Thus, such methods may result in mis- 
leading results. 
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FIG. 3: (color online) Comparison of the flow pattern in 
the reaction plane in 1.38A + 1.38A TeV peripheral, Pb+Pb 
collisions at two impact parameters: b — 0.56 max (left col- 
umn) and b = 0.7& max (right column), at a late stage of 10.16 
fm/c (240 time-steps). By this time an instability wave is 
also noticable at b = 0.5& max . The resolution increases from 
the top to the bottom as dx = 0.585, 0.4375, 0.35 fm (i.e. 
Nr — Rph/dx = 12, 16, 20). The number of marker particles 
were 8 3 , 6 3 and 5 3 for these resolutions, so that each of the 
marker particles carried about the same amount of baryon 
charge: 5.61, 5.61, 4.97 xl0 _J . For the impact parameter 
b = 0.56 m ax independently of the resolution and numerical 
viscosity the rotation of the dividing plane is 2 degrees, and 
the amplitude of the KH instability wave is not bigger than 
0.35 fm. For the impact parameter b = 0.7b max the rotation 
of the dividing plane is 3.75 degrees for all resolutions, but 
the amplitude of the KH instability wave is increasing with 
increasing resolution (decreasing numerical viscosity) as: 0.7, 
0.9 and 1.1 fm. 



A. CFD Results 

In heavy ion reactions the radial expansion modifies 
the dynamics of the development of the KH instability, 
but in case of low viscosity this is a strong instability 
and its initial signs can be clearly recognized in CFD 
calculations if both the viscosity and the numerical vis- 
cosity are sufficiently small. Looking at Fig. [T] initially 
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=0, the fluctuation is not visible, although the 
randomly placed markers include a possibility for fluctu- 



0.50 

0.45 
0.40 
0.35 
0.30 
0.25 
0.20 
0.15 
0.10 
0.05 
0.00 



Velocity in z-direction (c) 




. 4 




; I 40 




// 80 




;/ / - — ■ — . IzU 




if / ,■ 160 




/ t / 


200 








240 


- r/ S S 






y=0, z=0 



1 2 3 4 5 6 
x(fm) 



7 8 9 10 



cyc- 



FIG. 4: (color online) The velocity profile in the beam direc- 
tion as a function of the a>coordinate, at different times, 4, 
40, 80, 120, 160, 200 and 240 time-steps of 0.04233 fm/c for 
b — 0.76 max , cell size dx — 0.35 fm and 7 3 marker particles 
per fluid cell. The velocity is plotted in the reaction plane, (at 
y — and z = 0), and it is antisymmetric in the ±x-direction, 
but only the upper x > half of the reaction plane is plotted. 
In the CFD calculation the mirror symmetry is exact. 



ations (the length of the system is 35 cells). At 1.5 fm/c 
(-^Vcyc=35), there appear two semi-sinusoidal waves. The 
length of the system is 45 cells. The amplitude of tur- 
bulence is w 1 cell. At 3.0 fm/c (A cyc =70) the length 
of the system is « 50 cells. The central perturbation 
of 15 cell wavelength weakening while the outside part 
grows. The amplitude of turbulence is w 3 cells. At 4.4 
fm/c (iV C y C =105) the middle middle perturbation wave is 
weakening further, and the amplitude of the turbulence 
is w 4 cells. By 5.9 - 7.4 fm/c (7V cyc =140 H- 175) the mid- 
dle, short wavelength perturbation is hardly visible. The 
amplitude of the turbulence has reached 6 -j- 10 cells. 

Due to the fact that the radial expansion and the shear 
flow are superimposed upon each other the growth rate 
and the wave shape of the developing turbulence are not 
identifiable fully as in (quasi-) stationary flow. The ini- 
tial shorter wavelength perturbations become unstable 
and disappear, while the longest one grows due to the 
radial expansion. This can be attributed to the fact that 
the Kolmogorov minimal size is increasing faster than 
the expansion, thus the short perturbations are becom- 
ing unstable while the largest one survives. 

The phenomenon enables us to draw some quantita- 
tive consequences from the physical viscosity of QGP on 
qualitative differences in the flow pattern. The smaller 
central perturbation which disappears during expansion 
is not detectable. The larger one develops if the avail- 
able length of the system exceeds the Kolmogorov length 
scale. Then the shortest of these possible perturbations 
will grow fastest and will lead to enhanced and observable 
"rotation" . Increased beam energy leads to increasing V, 
which leads to an exponential increase in the growth of 
the instability, much more than just the linearly increased 
angular momentum would cause! 
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As the estimates of the previous section indicate, the 
development of the KH instability is critically dependent 
on the flow configuration, and just as the numerical es- 
timates for AkoI indicate the CFD results for impact pa- 
rameter b = 0.56 max also show that the KH instability is 
weak and it does not develop with increasing resolution 
for more central collisions, see Fig. [5] 

The dynamics in the CFD model indicates that the KH 
instability may start to develop in ultra-relativistic heavy 
ion collisions. The numerical viscosity in the calculation 
is about the same as the conjectured minimal viscosity 
@ as discussed in ref. [Icj . 

The analysis in the previous section assumed an initial 
viscosity based on ref. Q which was around the minimal 
viscosity and it was sufficient to develop the critical size 
KH instability. An order of magnitude larger viscosity 
would not be able to create a sufficiently small initial 
perturbation. 



B. Instability estimates for viscous fluids 

In shear flow one of the necessary conditions for the 
start a KH instability is that there should be an inflection 
point in the basic velocity profile u z (x) , or in other words 
u"(x) must change sign at least once @ Sect. 8. The 
initial condition |7J and the subsequent flow satisfy this 
requirement. 

As Fig. [4] shows, at time-step 4 and 40 the condition is 
satisfied in the center (x = 0), later at steps 80, 120 and 
160 a second wave develops, however the distance of the 
nodes is less than Ako1> so considering the viscous limits 
these secondary waves are not realizable physically. For a 
single wave of KH instability the situation is established, 
and persists up to time-step 120. By time-step 160 the 
central density drops considerably and apart from short 
length fluctuations the required condition is not satisfied. 

This is due to the fact that in the discussed heavy 
ion collision, we do not have a stationary boundary con- 
dition, and the spherical expansion competes with the 
growth of KH instability. At later stages, time-step 160 
and later, the growth speed, k V decreases, as both the 
velocity and the wave number decrease, while the expan- 
sion speed is increasing. Thus at time-step 160 the flow 
goes over to a combined approximately radial expansion 
and rotation. On the other hand by this time the system 
is close to freeze out and hadronization. At this time the 
matter is dilute and weakly interacting, so the local fluid 
dynamical equilibrium cannot be maintained and the FD 
description will not be applicable. 

The central zones at late stages have low density and 
low pressure, which resembles a cavitation bubble in the 
QGP liquid and this bubble will contain the condensed 
hadrons. Due to the expansion the bubble will not re- 
collapse as usual in cavitation phenomena in classical 
fluid flow; instead the surrounding QGP fluid will break 
up into peaces and will also hadronize and freeze out si- 
multaneously. 
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FIG. 5: (color online) The detailed view of the marker parti- 
cle positions in the lower half of the initial state markers after 
175 time-steps. A 1.38 A + 1.38 A TeV energy Pb+Pb periph- 
eral collision is shown, at b — 0.7 b max impact parameter with 
7 3 = 343 markers per initial, normal density fluid cell resolu- 
tion. The lines across the collision center point indicate the 
initial dividing axis, the change of this axis due to rotation 
and the additional change of rotation arising from the start- 
up of a Kelvin-Helmholtz type of instability. This additional 
effect more than doubles the rotation. In this calculation the 
cell size is dx — dy — dz — 0.35 fm, with a total number of 
1814814 marker particles. 



IV. CONCLUSIONS 

Based on theoretical estimates and on CFD model 
calculations, one should explore the possibility at LHC 
energies of a KH instability developing in peripheral, 
b = 0.6 — 0.86 max Pb+Pb collisions. The formation of 
the instability may take place up to about 5 fm/c, beyond 
which the radial expansion becomes dominant although 
the system still rotates. 

The KH instability is rather sensitive to the value of 
viscosity, so it is a perfect tool to measure the viscosity of 
QGP. The rotation of the weak anti-flow peak to forward 
angles was predicted earlier @. The rotation of the Vi- 
peak to forward angles depends sensitively on the balance 
between the speed of radial expansion/explosion and the 
initial angular momentum (which increases with increas- 
ing beam energy). If the radial expansion is stronger than 
estimated @ then the peak may remain an "antiflow" 
peak and the KH instability would destructively inter- 
fere with this peak. If the peak has rotated to forward 
angles as predicted in @ and used also in these calcu- 
lations, the KH instability increases the rotation and it 
converts a larger part of beam energy into rotation than 
it would happen in a simple solid body type of rotation. 
See Fig. El 

At cycle 160 (t = 6.77 fm/c) using the method of ref. 
Q primary v\ (y) values for massless pions were evaluated 
for two different impact parameters, b = 0.5, 0.76 max , and 
for different grid resolutions. The number of marker par- 
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tides were chosen so that the number of marker particles 
per baryon charge was about the same, 20120 for cell size 
dx = 0.35 fm and 17825 for cell size dx = 0.585 fm. For 
the Pb+Pb reaction at b = 0.56 max and dx — 0.585 fm 
resolution the directed flow peak was at the rapidity bin, 
y = 0.45 with a peak value of ui(0.45) = 0.177. This was 
taken to be 100%. By decreasing the cell size to dx — 0.35 
fm the peak value increased to ui(0.45) = 0.200, i.e. by 
13%. The changes of the directed flow peak values are 
shown in Table I. 

TABLE I. Change of Directed Flow 



Pb+Pb Pb+Pb Pb+Pb Pb+Pb 





b/b max 


0.5 


0.5 


0.7 


0.7 




dx,dy, dz [fm] 


0.585 


0.350 


0.585 


0.350 


0.35 


My) {%) 


-3 


10 


25 


47 


0.45 


vi(y) [%} 




13 


32 


53 


0.55 


My) [%} 


-3 


9 


20 


50 



For b = 0.56 max the KH instability on the vi-peak is 
weak, 13%, which can partly be attributed to the de- 
creased viscous dissipation. At b = 0.7& max , the increase 
is 21%, it is significantly stronger. These primary data of 
of course are reduced by random initial state dissipation, 
as discussed in ref. [6j. The position of the peak in ra- 
pidity is hardly changing with higher grid resolution, for 
b = 0.56 max it moves from y — 0.485 to y = 0.46, while 
for b = 0.7& max it moves from y = 0.475 to y = 0.46. 

Recently if was pointed out 15] that due to random ini- 
tial fluctuations turbulence may show up and even grow 



in the transverse, [x, y] plane. The energy of the growth 
is provided by absorbing small, higher wave number per- 
turbations. We also observed this effect (see Fig. [1]). In 
that work an alternative detection method is suggested 
via measuring two particle correlations, which may also 
be used to detect the KH instability in the reaction plane. 

Although, the predicted rotation effect is not easily 
detectable due to initial state fluctuations, the KH in- 
stability enhances the flow and changes its pattern in 
peripheral collisions. The present developments suggest 
that the global collective v\ flow can be disentangled from 
random fluctuations. This is necessary to measure the 
global collective flow in peripheral collisions. The oppo- 
site, the separation of the flow originating from random 
initial state fluctuations is done successfully recently [l6[ 
for selected central collisions. 

The KH instability is very sensitive to the magnitude 
of the viscosity. Thus if this research is successful the 
analysis of global collective v\ flow as a function of beam 
energy and impact parameter may provide a precision 
measurement of viscosity and its variation. 
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